load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"


;===============================================================================
; Begin of Main Process
;===============================================================================

;;; read wind measurement data ;;;

csvFile = "/public/users/wrf/Projects/xiangshui_rttower.csv"

lines = readAsciiTable(csvFile,1,"string",1)
delim = ","

towerid_m = strip_quotes ( ndtooned (str_get_field(lines, 1, delim)))
timeStr_m   = strip_quotes ( ndtooned ( str_get_field (lines,2,delim) ) )
windsp10_m  = tofloat ( strip_quotes ( ndtooned ( str_get_field (lines,3,delim) ) ) )
winddir10_m = tointeger ( strip_quotes ( ndtooned ( str_get_field (lines,4,delim) ) ) )
windsp70_m  = tofloat ( strip_quotes ( ndtooned ( str_get_field (lines,33,delim) ) ) )
winddir70_m = tointeger ( strip_quotes ( ndtooned ( str_get_field (lines,34,delim) ) ) )

time_udunit_m = timestr_to_udunits(timeStr_m)
timeStr2_m = ut_string(time_udunit_m, "%Y-%N-%D %H:%M:%S")

; find the first and last day
indexLastMeasure = ind_nearest_coord( timestr_to_udunits("2012-11-08 00:00:00"), time_udunit_m, 0)
print(indexLastMeasure)
print(tostring(indexLastMeasure) + " " + tostring(time_udunit_m(indexLastMeasure)) + " " + timeStr2_m(indexLastMeasure))
print(timeStr2_m(0))

; save result to csv file
firstTimeRec = time_udunit_m(0)
lastTimeRec  = time_udunit_m(indexLastMeasure)
validTimes   = fspan(firstTimeRec, lastTimeRec, tointeger((lastTimeRec - firstTimeRec) * 2 + 1))
validIndex   = ind_nearest_coord(validTimes, time_udunit_m, 0)

validTowerid  = towerid_m(validIndex)
validTimeStrs = timeStr2_m(validIndex)
validWindsp10 = windsp10_m(validIndex)
validWinddir10 = winddir10_m(validIndex)
validWindsp70 = windsp70_m(validIndex)
validWinddir70 = winddir70_m(validIndex)
printArray = ( validTowerid + ", " + validTimeStrs \
                            + ", "+ sprintf("%6.2f", validWindsp10) + ", "+ sprinti("%3i", validWinddir10) \
                            + ", "+ sprintf("%6.2f", validWindsp70) + ", "+ sprinti("%3i", validWinddir70))

; write to csv file
output_meas_csvfile = "~/Projects/xiangshui_meas.csv"
asciiwrite(output_meas_csvfile, printArray)

delete(printArray)

;;; read NWP predicted wind data ;;;

; get nwp output txt file list, prediction time start at 20:00
nwpTxtPath        = "/public/users/wrf/wrfpost/test/station_wind"
nwpTxtFileList_20 = systemfunc("ls " + nwpTxtPath + "/wind_output_cixi_d02_*_12.txt")
nwpTxtFileList_08 = systemfunc("ls " + nwpTxtPath + "/wind_output_cixi_d02_*_00.txt")
nFile             = dimsizes(nwpTxtFileList_20)

delayTime_20 = 28
delayTime_08 = 16
dt           = 0.5 ; time resolution: 30min
allRecSize   = tointeger(nFile*24/dt)
timeStr_p    = new(allRecSize, "string")
windsp10_p   = new(allRecSize, "float")
winddir10_p  = new(allRecSize, "integer")
windsp70_p   = new(allRecSize, "float")
winddir70_p  = new(allRecSize, "integer")

; read wind data from each file
do ifile = 0, dimsizes(nwpTxtFileList_08) - 1
    delete(lines)
    lines = readAsciiTable(nwpTxtFileList_08(ifile),1,"string",1)

    timeStr_t   = ndtooned( str_get_field(lines, 3, delim))
    windsp10_t  = ndtooned( tofloat( str_get_field(lines,4,delim)))
    winddir10_t = ndtooned( tointeger( str_get_field(lines,5,delim)))
    windsp70_t  = ndtooned( tofloat( str_get_field(lines,6,delim)))
    winddir70_t = ndtooned( tointeger( str_get_field(lines,7,delim)))

    ; convert time string to udunit
    time_udunit_t = timestr_to_udunits(timeStr_t)

    ; use time as coordinate of variables
    windsp10_t!0     = "Time"
    windsp70_t!0     = "Time"
    winddir10_t!0    = "Time"
    winddir70_t!0    = "Time"
    windsp10_t&Time  = time_udunit_t
    windsp70_t&Time  = time_udunit_t
    winddir10_t&Time = time_udunit_t
    winddir70_t&Time = time_udunit_t

    ; find the period of 24h prediction
    predictionStartTime = time_udunit_t(0)
    firstValidTime      = predictionStartTime + delayTime_08
    lastValidTime       = predictionStartTime + delayTime_08 + 24 - dt
    copy_VarAtts(predictionStartTime, firstValidTime)
    copy_VarAtts(predictionStartTime, lastValidTime)
    selectTimes         = fspan(firstValidTime,lastValidTime,24*2)
    nRecDaily           = dimsizes(selectTimes)
    index_selectedTime  = ind_nearest_coord(selectTimes, time_udunit_t, 0)
    windsp10_p(ifile*nRecDaily:((ifile+1)*nRecDaily-1)) = windsp10_t(index_selectedTime)
    windsp70_p(ifile*nRecDaily:((ifile+1)*nRecDaily-1)) = windsp70_t(index_selectedTime)
    winddir10_p(ifile*nRecDaily:((ifile+1)*nRecDaily-1)) = winddir10_t(index_selectedTime)
    winddir70_p(ifile*nRecDaily:((ifile+1)*nRecDaily-1)) = winddir70_t(index_selectedTime)
    timeStr_p(ifile*nRecDaily:((ifile+1)*nRecDaily-1)) = timeStr_t(index_selectedTime)

    delete(timeStr_t)
    delete(time_udunit_t)
    delete(windsp10_t)
    delete(winddir10_t)
    delete(windsp70_t)
    delete(winddir70_t)
end do

printArray = ("41101," +timeStr_p \
                    + ", " + sprintf("%6.2f", windsp10_p) +", "+ sprinti("%3i", winddir10_p) \
                    + ", " + sprintf("%6.2f", windsp70_p) +", "+ sprinti("%3i", winddir70_p))

; write to csv file
output_meas_csvfile = "~/Projects/cixi_predict_24_08.csv"
asciiwrite(output_meas_csvfile, printArray)
